clear
clc

%choose a competition
zz=5;

%load data
w=load('competitions.mat');
COMP=w.comp_list(zz,1);
x=load('data.mat');
y=load('qestimates.mat');
competitionid=x.competitionid;
maxscore = x.maxscore(competitionid==COMP,1);
pubscore_normal = x.pubscore_normal(competitionid==COMP,1);
length_days= x.length_days(competitionid==COMP,1);
length_days=length_days(1);
rewardquantity= x.rewardquantity(competitionid==COMP,1);
rewardquantity=rewardquantity(1);
mergers_data= x.totmergers(competitionid==COMP,1);
mergers_data=mergers_data(1);
nSim = 500;

%%%%%state space parameters
N=40; %number of individuals
T=360; %number of periods
Svalues=unique([unique(maxscore);max(maxscore)+[0.002:0.002:0.08]']);
nScore = 20;
Svalues = Svalues(3:nScore+2,1);
nsubs_data=size(maxscore,1);
lambda1=0.75;
lambda2=1-lambda1;

%%%%%q-function
betachangehat=[y.beta0(y.cid==COMP,1) ; y.beta1(y.cid==COMP,1); y.beta1tr(y.cid==COMP,1)];
probchange = @(score) (exp(betachangehat(1) + betachangehat(2)*score)./(1+exp(betachangehat(1) + betachangehat(2)*score)));
q_sp=probchange(Svalues);
%q_sp(end,1)=0;
probchange = @(score) (exp(betachangehat(1) + betachangehat(3) + betachangehat(2)*score)./(1+exp(betachangehat(1)+ betachangehat(3) + betachangehat(2)*score)));
q_team=probchange(Svalues);
%q_team(end,1)=0;

%load estimates
load(sprintf('Estimates/%s_%02d.mat', 'estimates',COMP))

%compute equilibrium
%Values
L_sp_end=1;
L_team_end=1/2;
F_end = 0;
L_sp=zeros(N/2+1,nScore,T);
F_sp=zeros(N/2+1,nScore,T,2);
L_team=zeros(N/2+1,nScore,T);
F_team=F_sp;
psi=1-lambda1-lambda2+2*lambda2*[0:1:N/2]'/N;
Psub_f_team=F_sp;
Psub_f_sp=F_sp;
Psub_l_team=L_sp;
Psub_l_sp=L_sp;
Pteam=F_sp;

L_team(:,:,end)=L_team_end;
L_sp(:,:,end)=L_sp_end;


%%

K = @(x,sigma) x.^sigma;
EC_Z = @(z,sigma) max(z,0)*sigma/(1+sigma); %E[c|c<z] (truncated expectation)
fun_p = @(q,V1,V2,sigma) K(max(q*(V1-V2),0),sigma);
fun_exp_max = @(p,q,V1,V2,sigma) p*(q*V1+(1-q)*V2-EC_Z(q*(V1-V2),sigma))+(1-p)*V2;

% for k=1:T-1
%     t = T - k;
%     for i=0:(N/2)
%         aux_Teams = i;
%         iaux_Teams = aux_Teams+1;
%         aux_sp=N - 2*aux_Teams;
% %         L_sp(:,nScore,t)=L_sp_end;
% %         L_team(:,nScore,t)=L_team_end;
%         for j=1:nScore-1
%             aux_Score = j;
%             
%             p_team_f_sp_0=0;
%             p_team_f_sp_1=0;
%             F_sp_forms_team_0=0;
%             F_sp_forms_team_1=0;
%             F_sp_rival_sp_0=0;
%             F_sp_rival_sp_1=0;
%             F_sp_rival_team_0=0;
%             F_sp_rival_team_1=0;
%             F_sp_team_forms_0=0;
%             F_sp_team_forms_1=0;
%             L_sp_team_forms=0;
%             L_team_team_forms=0;
%             F_team_rival_team_0=0;
%             F_team_rival_team_1=0;
%             F_team_rival_sp_0=0;
%             F_team_rival_sp_1=0;
%             F_team_team_forms_0=0;
%             F_team_team_forms_1=0;
%             
%              %%%%%%%%CCPs
%             p_sub_f_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaS);
%             p_sub_f_sp_1 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaS);
%             p_sub_f_team_0 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
%             p_sub_f_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
%             p_sub_l_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
%             p_sub_l_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
%             if aux_Teams<N/2
%                 p_team_f_sp_0 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaT);
%                 p_team_f_sp_1 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,2),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaT);
%                 
%              %aux_F_sp_forms_team computes F_sp_forms team as a function of
%             %F^{sp}_{s,n,l}, l, and p^team
%                 aux_F_sp_forms_team = @(p,l) fun_exp_max(p,1,F_team(iaux_Teams+1,aux_Score,t+1,1+l),F_sp(iaux_Teams,aux_Score,t+1,1+l),sigmaT);
%                 F_sp_forms_team_0 = aux_F_sp_forms_team(p_team_f_sp_0,0);
%                 F_sp_forms_team_1 = aux_F_sp_forms_team(p_team_f_sp_1,1);                
%  
%                 %computing F^{sp,rival sp}
%                 aux_F_sp_rival_sp= @(pl,pf,F,l) ((1-l)/(aux_sp-1))*(pl*(q_sp(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-(1-l))/(aux_sp-1))*(pf*(q_sp(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F)+(1-pf)*F);
%                 F_sp_rival_sp_0 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
%                 F_sp_rival_sp_1 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
%                 L_sp_team_forms = p_team_f_sp_0*L_sp(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
%                 L_team_team_forms = p_team_f_sp_1*L_team(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
%                
%                 aux_F_team_rival_sp= @(pl,pf,l) ((1-l)/(aux_sp))*(pl*(q_sp(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-1+l)/(aux_sp))*(pf*(q_sp(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
%                 F_team_rival_sp_0 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,0);
%                 F_team_rival_sp_1 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,1);
%                 F_team_team_forms_0 = p_team_f_sp_0*F_team(iaux_Teams+1,aux_Score,t+1,1+0) + (1-p_team_f_sp_0)*F_team(iaux_Teams,aux_Score,t+1,1+0);
%                 F_team_team_forms_1 = p_team_f_sp_1*F_team(iaux_Teams+1,aux_Score,t+1,1+1) + (1-p_team_f_sp_1)*F_team(iaux_Teams,aux_Score,t+1,1+1);
%                 
%             end
%             
%               
%                 
%             %%%%%%%%%F_sp
%             %aux_F_sp_own computes F_sp_own as a function of F^{sp}_{s,n,l} and
%             %p^sub
%             aux_F_sp_own = @(p,F) fun_exp_max(p,q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),F,sigmaS);
%             F_sp_own_0 = aux_F_sp_own(p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1));
%             F_sp_own_1 = aux_F_sp_own(p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,2));
%             
% 
%             %computing F^{sp,rival team}
%             if aux_Teams>0
%                 aux_F_sp_rival_team = @(pl,pf,F,l) (l/aux_Teams)*(pl*(q_team(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l)/aux_Teams)*(pf*(q_team(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F)+(1-pf)*F);
%                 F_sp_rival_team_0 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
%                 F_sp_rival_team_1 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
% 
%             end
%             
%             %Computing F_sp_team forms
%             if aux_sp>2
%                 aux_F_sp_team_forms = @(p,l) ((1-p)*((aux_sp-2+l)/(aux_sp-1))+((1-l)/(aux_sp-1)))*F_sp(iaux_Teams,aux_Score,t+1,1+l) + p*((aux_sp-2+l)/(aux_sp-1))*((1/(aux_sp-2+l))*F_team(iaux_Teams+1,aux_Score,t+1,1+l) + ((aux_sp-3+l)/(aux_sp-2+l))*F_sp(iaux_Teams+1,aux_Score,t+1,1+l));
%                 F_sp_team_forms_0 = aux_F_sp_team_forms(p_team_f_sp_0,0);
%                 F_sp_team_forms_1 = aux_F_sp_team_forms(p_team_f_sp_1,1);
%             end
%             
%             F_sp(iaux_Teams,aux_Score,t,1) = (lambda1/N)*F_sp_own_0 + (lambda2/N)*F_sp_forms_team_0 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,1) + (lambda1*(2*aux_Teams)/N)*F_sp_rival_team_0 + lambda1*((aux_sp-1)/N)*F_sp_rival_sp_0 + lambda2*((aux_sp-1)/N)*F_sp_team_forms_0;
%             F_sp(iaux_Teams,aux_Score,t,2) = (lambda1/N)*F_sp_own_1 + (lambda2/N)*F_sp_forms_team_1 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,2) + (lambda1*(2*aux_Teams)/N)*F_sp_rival_team_1 + lambda1*((aux_sp-1)/N)*F_sp_rival_sp_1 + lambda2*((aux_sp-1)/N)*F_sp_team_forms_1;
%             
%             %%%%%%%%%L_sp
%             L_sp_own_0 = fun_exp_max(p_sub_l_sp_0,q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
%             L_sp_rival_team = p_sub_f_team_0*(q_team(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_0)*L_sp(iaux_Teams,aux_Score,t+1);
%             L_sp_rival_sp = p_sub_f_sp_0*(q_sp(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
% 
%             L_sp(iaux_Teams,aux_Score,t) = (lambda1/N)*L_sp_own_0 + (psi(iaux_Teams)+lambda2/N)*L_sp(iaux_Teams,aux_Score,t+1) + (lambda1*(2*aux_Teams)/N)*L_sp_rival_team + lambda1*((aux_sp-1)/N)*L_sp_rival_sp + lambda2*((aux_sp-1)/N)*L_sp_team_forms;
%             
%             
%             %%%%%%%%%L_team
%             L_team_own_1 = fun_exp_max(p_sub_l_team_1,q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
%             L_team_rival_team = p_sub_f_team_1*(q_team(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_1)*L_team(iaux_Teams,aux_Score,t+1);
%             L_team_rival_sp = p_sub_f_sp_1*(q_sp(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
% 
%             L_team(iaux_Teams,aux_Score,t) = ((2*lambda1)/N)*L_team_own_1 + psi(iaux_Teams)*L_team(iaux_Teams,aux_Score,t+1) + (lambda1*(2*(aux_Teams-1))/N)*L_team_rival_team + lambda1*(aux_sp/N)*L_team_rival_sp + lambda2*(aux_sp/N)*L_team_team_forms;
%             
%             
%             %%%%%%%%%%F_team
%             f_team_own_0 = fun_exp_max(p_sub_f_team_0,q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
%             f_team_own_1 = fun_exp_max(p_sub_f_team_1,q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
%             
%             %computing F^{team,rival team}
%             if aux_Teams>1
%                 aux_F_team_rival_team = @(pl,pf,l) (l/(aux_Teams-1))*(pl*(q_team(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l-1)/(aux_Teams-1))*(pf*(q_team(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
%                 F_team_rival_team_0 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_0,0);
%                 F_team_rival_team_1 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_1,1);
%             end
%             
%     
%             F_team(iaux_Teams,aux_Score,t,1+0) = ((2*lambda1)/N)*f_team_own_0 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+0) + (lambda1*(2*(aux_Teams-1))/N)*F_team_rival_team_0 + lambda1*((aux_sp)/N)*F_team_rival_sp_0 + lambda2*((aux_sp)/N)*F_team_team_forms_0;
%             F_team(iaux_Teams,aux_Score,t,1+1) = ((2*lambda1)/N)*f_team_own_1 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+1) + (lambda1*(2*(aux_Teams-1))/N)*F_team_rival_team_1 + lambda1*((aux_sp)/N)*F_team_rival_sp_1 + lambda2*((aux_sp)/N)*F_team_team_forms_1;
%                        
%             %%%% Storing CCPs
%             Psub_f_team(iaux_Teams,aux_Score,t,1)=p_sub_f_team_0;
%             Psub_f_team(iaux_Teams,aux_Score,t,2)=p_sub_f_team_1;
%             Psub_f_sp(iaux_Teams,aux_Score,t,1)=p_sub_f_sp_0;
%             Psub_f_sp(iaux_Teams,aux_Score,t,2)=p_sub_f_sp_1;
%             Psub_l_team(iaux_Teams,aux_Score,t)=p_sub_l_team_1;
%             Psub_l_sp(iaux_Teams,aux_Score,t)=p_sub_l_sp_0;
%             Pteam(iaux_Teams,aux_Score,t,1)=p_team_f_sp_0;
%             Pteam(iaux_Teams,aux_Score,t,2)=p_team_f_sp_1;
%                    
%         end
%     end
% [k T];    
% end

% K = @(x,sigma) x.^sigma;
% EC_Z = @(z,sigma) z*sigma/(1+sigma); %E[c|c<z] (truncated expectation)
% fun_p = @(q,V1,V2,sigma) K(max(q*(V1-V2),0),sigma);
% fun_exp_max = @(p,q,V1,V2,sigma) p*(q*V1+(1-q)*V2-EC_Z(q*(V1-V2),sigma))+(1-p)*V2;
% 
% for k=1:T-1
%     t = T - k;
%     for i=0:(N/2)
%         aux_Teams = i;
%         iaux_Teams = aux_Teams+1;
%         aux_sp=N - 2*aux_Teams;
% %         L_sp(:,nScore,t)=L_sp_end;
% %         L_team(:,nScore,t)=L_team_end;
%         for j=1:nScore
%             aux_Score = j;
%             
%             p_team_f_sp_0=0;
%             p_team_f_sp_1=0;
%             F_sp_forms_team_0=0;
%             F_sp_forms_team_1=0;
%             F_sp_rival_sp_0=0;
%             F_sp_rival_sp_1=0;
%             F_sp_rival_team_0=0;
%             F_sp_rival_team_1=0;
%             F_sp_team_forms_0=0;
%             F_sp_team_forms_1=0;
%             L_sp_team_forms=0;
%             L_team_team_forms=0;
%             F_team_rival_team_0=0;
%             F_team_rival_team_1=0;
%             F_team_rival_sp_0=0;
%             F_team_rival_sp_1=0;
%             F_team_team_forms_0=0;
%             F_team_team_forms_1=0;
%             
%              %%%%%%%%CCPs
%             p_sub_f_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaS);
%             p_sub_f_sp_1 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaS);
%             p_sub_f_team_0 = fun_p(q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
%             p_sub_f_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
%             p_sub_l_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
%             p_sub_l_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
%             if aux_Teams<N/2
%                 p_team_f_sp_0 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaT);
%                 p_team_f_sp_1 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,2),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaT);
%                 
%              %aux_F_sp_forms_team computes F_sp_forms team as a function of
%             %F^{sp}_{s,n,l}, l, and p^team
%                 aux_F_sp_forms_team = @(p,l) fun_exp_max(p,1,F_team(iaux_Teams+1,aux_Score,t+1,1+l),F_sp(iaux_Teams,aux_Score,t+1,1+l),sigmaT);
%                 F_sp_forms_team_0 = aux_F_sp_forms_team(p_team_f_sp_0,0);
%                 F_sp_forms_team_1 = aux_F_sp_forms_team(p_team_f_sp_1,1);                
%  
%                 %computing F^{sp,rival sp}
%                 aux_F_sp_rival_sp= @(pl,pf,F,l) ((1-l)/(aux_sp-1))*(pl*(q_sp(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-(1-l))/(aux_sp-1))*(pf*(q_sp(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F)+(1-pf)*F);
%                 F_sp_rival_sp_0 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
%                 F_sp_rival_sp_1 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
%                 L_sp_team_forms = p_team_f_sp_0*L_sp(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
%                 L_team_team_forms = p_team_f_sp_1*L_team(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
%                
%                 aux_F_team_rival_sp= @(pl,pf,l) ((1-l)/(aux_sp))*(pl*(q_sp(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-1+l)/(aux_sp))*(pf*(q_sp(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
%                 F_team_rival_sp_0 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,0);
%                 F_team_rival_sp_1 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,1);
%                 F_team_team_forms_0 = p_team_f_sp_0*F_team(iaux_Teams+1,aux_Score,t+1,1+0) + (1-p_team_f_sp_0)*F_team(iaux_Teams,aux_Score,t+1,1+0);
%                 F_team_team_forms_1 = p_team_f_sp_1*F_team(iaux_Teams+1,aux_Score,t+1,1+1) + (1-p_team_f_sp_1)*F_team(iaux_Teams,aux_Score,t+1,1+1);
%                 
%             end
%             
%               
%                 
%             %%%%%%%%%F_sp
%             %aux_F_sp_own computes F_sp_own as a function of F^{sp}_{s,n,l} and
%             %p^sub
%             aux_F_sp_own = @(p,F) fun_exp_max(p,q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),F,sigmaS);
%             F_sp_own_0 = aux_F_sp_own(p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1));
%             F_sp_own_1 = aux_F_sp_own(p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,2));
%             
% 
%             %computing F^{sp,rival team}
%             if aux_Teams>0
%                 aux_F_sp_rival_team = @(pl,pf,F,l) (l/aux_Teams)*(pl*(q_team(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l)/aux_Teams)*(pf*(q_team(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F)+(1-pf)*F);
%                 F_sp_rival_team_0 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
%                 F_sp_rival_team_1 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
% 
%             end
%             
%             %Computing F_sp_team forms
%             if aux_sp>2
%                 aux_F_sp_team_forms = @(p,l) ((1-p)*((aux_sp-2+l)/(aux_sp-1))+((1-l)/(aux_sp-1)))*F_sp(iaux_Teams,aux_Score,t+1,1+l) + p*((aux_sp-2+l)/(aux_sp-1))*((1/(aux_sp-2+l))*F_team(iaux_Teams+1,aux_Score,t+1,1+l) + ((aux_sp-3+l)/(aux_sp-2+l))*F_sp(iaux_Teams+1,aux_Score,t+1,1+l));
%                 F_sp_team_forms_0 = aux_F_sp_team_forms(p_team_f_sp_0,0);
%                 F_sp_team_forms_1 = aux_F_sp_team_forms(p_team_f_sp_1,1);
%             end
%             
%             F_sp(iaux_Teams,aux_Score,t,1) = (lambda1/N)*F_sp_own_0 + (lambda2/N)*F_sp_forms_team_0 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,1) + (lambda1*(2*aux_Teams)/N)*F_sp_rival_team_0 + lambda1*((aux_sp-1)/N)*F_sp_rival_sp_0 + lambda2*((aux_sp-1)/N)*F_sp_team_forms_0;
%             F_sp(iaux_Teams,aux_Score,t,2) = (lambda1/N)*F_sp_own_1 + (lambda2/N)*F_sp_forms_team_1 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,2) + (lambda1*(2*aux_Teams)/N)*F_sp_rival_team_1 + lambda1*((aux_sp-1)/N)*F_sp_rival_sp_1 + lambda2*((aux_sp-1)/N)*F_sp_team_forms_1;
%             
%             %%%%%%%%%L_sp
%             L_sp_own_0 = fun_exp_max(p_sub_l_sp_0,q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
%             L_sp_rival_team = p_sub_f_team_0*(q_team(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_0)*L_sp(iaux_Teams,aux_Score,t+1);
%             L_sp_rival_sp = p_sub_f_sp_0*(q_sp(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
% 
%             L_sp(iaux_Teams,aux_Score,t) = (lambda1/N)*L_sp_own_0 + (psi(iaux_Teams)+lambda2/N)*L_sp(iaux_Teams,aux_Score,t+1) + (lambda1*(2*aux_Teams)/N)*L_sp_rival_team + lambda1*((aux_sp-1)/N)*L_sp_rival_sp + lambda2*((aux_sp-1)/N)*L_sp_team_forms;
%             
%             
%             %%%%%%%%%L_team
%             L_team_own_1 = fun_exp_max(p_sub_l_team_1,q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
%             L_team_rival_team = p_sub_f_team_1*(q_team(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_1)*L_team(iaux_Teams,aux_Score,t+1);
%             L_team_rival_sp = p_sub_f_sp_1*(q_sp(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
% 
%             L_team(iaux_Teams,aux_Score,t) = ((2*lambda1)/N)*L_team_own_1 + psi(iaux_Teams)*L_team(iaux_Teams,aux_Score,t+1) + (lambda1*(2*(aux_Teams-1))/N)*L_team_rival_team + lambda1*(aux_sp/N)*L_team_rival_sp + lambda2*(aux_sp/N)*L_team_team_forms;
%             
%             
%             %%%%%%%%%%F_team
%             f_team_own_0 = fun_exp_max(p_sub_f_team_0,q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
%             f_team_own_1 = fun_exp_max(p_sub_f_team_1,q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
%             
%             %computing F^{team,rival team}
%             if aux_Teams>1
%                 aux_F_team_rival_team = @(pl,pf,l) (l/(aux_Teams-1))*(pl*(q_team(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l-1)/(aux_Teams-1))*(pf*(q_team(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
%                 F_team_rival_team_0 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_0,0);
%                 F_team_rival_team_1 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_1,1);
%             end
%             
%     
%             F_team(iaux_Teams,aux_Score,t,1+0) = ((2*lambda1)/N)*f_team_own_0 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+0) + (lambda1*(2*(aux_Teams-1))/N)*F_team_rival_team_0 + lambda1*((aux_sp)/N)*F_team_rival_sp_0 + lambda2*((aux_sp)/N)*F_team_team_forms_0;
%             F_team(iaux_Teams,aux_Score,t,1+1) = ((2*lambda1)/N)*f_team_own_1 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+1) + (lambda1*(2*(aux_Teams-1))/N)*F_team_rival_team_1 + lambda1*((aux_sp)/N)*F_team_rival_sp_1 + lambda2*((aux_sp)/N)*F_team_team_forms_1;
%                        
%             %%%% Storing CCPs
%             Psub_f_team(iaux_Teams,aux_Score,t,1)=p_sub_f_team_0;
%             Psub_f_team(iaux_Teams,aux_Score,t,2)=p_sub_f_team_1;
%             Psub_f_sp(iaux_Teams,aux_Score,t,1)=p_sub_f_sp_0;
%             Psub_f_sp(iaux_Teams,aux_Score,t,2)=p_sub_f_sp_1;
%             Psub_l_team(iaux_Teams,aux_Score,t)=p_sub_l_team_1;
%             Psub_l_sp(iaux_Teams,aux_Score,t)=p_sub_l_sp_0;
%             Pteam(iaux_Teams,aux_Score,t,1)=p_team_f_sp_0;
%             Pteam(iaux_Teams,aux_Score,t,2)=p_team_f_sp_1;
%                    
%         end
%     end
% [k T];    
% end

for k=1:T-1
    t = T - k;
    for i=0:(N/2)
        aux_Teams = i;
        iaux_Teams = aux_Teams+1;
        aux_sp=N - 2*aux_Teams;
%         L_sp(:,nScore,t)=L_sp_end;
%         L_team(:,nScore,t)=L_team_end;
        for j=1:nScore
            aux_Score = j;
            
            %%%%%%%%CCPs
            p_team_f_sp_0=0;
            p_team_f_sp_1=0;            
            p_sub_f_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaS);
            p_sub_f_sp_1 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaS);
            p_sub_f_team_0 = fun_p(q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
            p_sub_f_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
            p_sub_l_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
            p_sub_l_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
            if aux_Teams<N/2
                p_team_f_sp_0 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaT);
                p_team_f_sp_1 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,2),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaT);
            end
            
            
            %%%%%%%%%F_sp
            %aux_F_sp_own computes F_sp_own as a function of F^{sp}_{s,n,l} and
            %p^sub
            aux_F_sp_own = @(p,F) fun_exp_max(p,q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),F,sigmaS);
            F_sp_own_0 = aux_F_sp_own(p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1));
            F_sp_own_1 = aux_F_sp_own(p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,2));
            
            %aux_F_sp_forms_team computes F_sp_forms team as a function of
            %F^{sp}_{s,n,l}, l, and p^team
            if aux_Teams<N/2
                aux_F_sp_forms_team = @(p,l) fun_exp_max(p,1,F_team(iaux_Teams+1,aux_Score,t+1,1+l),F_sp(iaux_Teams,aux_Score,t+1,1+l),sigmaT);
                F_sp_forms_team_0 = aux_F_sp_forms_team(p_team_f_sp_0,0);
                F_sp_forms_team_1 = aux_F_sp_forms_team(p_team_f_sp_1,1);
            else
                F_sp_forms_team_0=0;
                F_sp_forms_team_1=0;
            end
            
            %computing F^{sp,rival sp}
            if aux_Teams<N/2
                aux_F_sp_rival_sp= @(pl,pf,F,l) ((1-l)/(aux_sp-1))*(pl*(q_sp(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-(1-l)-1)/(aux_sp-1))*(pf*(q_sp(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F)+(1-pf)*F);
                F_sp_rival_sp_0 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
                F_sp_rival_sp_1 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
            else
                F_sp_rival_sp_0=0;
                F_sp_rival_sp_1=0;
            end
            
            %computing F^{sp,rival team}
            if aux_Teams>0
                aux_F_sp_rival_team = @(pl,pf,F,l) (l/aux_Teams)*(pl*(q_team(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l)/aux_Teams)*(pf*(q_team(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F)+(1-pf)*F);
                F_sp_rival_team_0 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
                F_sp_rival_team_1 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
            else
                F_sp_rival_team_0=0;
                F_sp_rival_team_1=0;
            end
            
            %Computing F_sp_team forms
            if aux_sp>2
                aux_F_sp_team_forms = @(p,l) ((1-p)*((aux_sp-2+l)/(aux_sp-1))+((1-l)/(aux_sp-1)))*F_sp(iaux_Teams,aux_Score,t+1,1+l) + p*((aux_sp-2+l)/(aux_sp-1))*((1/(aux_sp-2+l))*F_team(iaux_Teams+1,aux_Score,t+1,1+l) + ((aux_sp-3+l)/(aux_sp-2+l))*F_sp(iaux_Teams+1,aux_Score,t+1,1+l));
                F_sp_team_forms_0 = aux_F_sp_team_forms(p_team_f_sp_0,0);
                F_sp_team_forms_1 = aux_F_sp_team_forms(p_team_f_sp_1,1);
            else
                F_sp_team_forms_0=0;
                F_sp_team_forms_1=0;
            end
            
            F_sp(iaux_Teams,aux_Score,t,1) = (lambda1/N)*F_sp_own_0 + (lambda2/N)*F_sp_forms_team_0 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,1) + (lambda1*(2*aux_Teams)/N)*max(0,F_sp_rival_team_0) + max(0,lambda1*((aux_sp-1)/N))*max(0,F_sp_rival_sp_0) + max(0,lambda2*((aux_sp-1)/N))*max(0,F_sp_team_forms_0);
            F_sp(iaux_Teams,aux_Score,t,2) = (lambda1/N)*F_sp_own_1 + (lambda2/N)*F_sp_forms_team_1 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,2) + (lambda1*(2*aux_Teams)/N)*max(0,F_sp_rival_team_1) + max(0,lambda1*((aux_sp-1)/N))*max(0,F_sp_rival_sp_1) + max(0,lambda2*((aux_sp-1)/N))*max(0,F_sp_team_forms_1);
            if aux_Teams==N/2
              F_sp(iaux_Teams,aux_Score,t,1)=0;
              F_sp(iaux_Teams,aux_Score,t,2)=0;
            end
            
            %%%%%%%%%L_sp
            L_sp_own_0 = fun_exp_max(p_sub_l_sp_0,q_sp(aux_Score),L_sp(iaux_Teams,min(aux_Score+1,nScore),t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
            L_sp_rival_team = p_sub_f_team_0*(q_team(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_0)*L_sp(iaux_Teams,aux_Score,t+1);
            L_sp_rival_sp = p_sub_f_sp_0*(q_sp(aux_Score)*F_sp(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
            if aux_Teams<N/2
                L_sp_team_forms = p_team_f_sp_0*L_sp(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
            else
                L_sp_team_forms=0;
            end
            L_sp(iaux_Teams,aux_Score,t) = (lambda1/N)*L_sp_own_0 + (psi(iaux_Teams)+lambda2/N)*L_sp(iaux_Teams,aux_Score,t+1) + (lambda1*(2*aux_Teams)/N)*max(0,L_sp_rival_team) + max(0,lambda1*((aux_sp-1)/N))*max(0,L_sp_rival_sp) + max(0,lambda2*((aux_sp-1)/N))*max(0,L_sp_team_forms);
            if aux_Teams==N/2
              L_sp(iaux_Teams,aux_Score,t)=0;
              L_sp(iaux_Teams,aux_Score,t)=0;
            end            
            
            %%%%%%%%%L_team
            L_team_own_1 = fun_exp_max(p_sub_l_team_1,q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
            L_team_rival_team = p_sub_f_team_1*(q_team(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_1)*L_team(iaux_Teams,aux_Score,t+1);
            L_team_rival_sp = p_sub_f_sp_1*(q_sp(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
            if aux_Teams<N/2
                L_team_team_forms = p_team_f_sp_1*L_team(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
            else
                L_team_team_forms=0;
            end
            L_team(iaux_Teams,aux_Score,t) = ((2*lambda1)/N)*L_team_own_1 + psi(iaux_Teams)*L_team(iaux_Teams,aux_Score,t+1) + max(0,(lambda1*(2*(aux_Teams-1))/N))*max(0,L_team_rival_team) + lambda1*(aux_sp/N)*max(0,L_team_rival_sp) + lambda2*(aux_sp/N)*max(0,L_team_team_forms);
            if aux_Teams==0
            L_team(iaux_Teams,aux_Score,t) = psi(iaux_Teams)*L_team(iaux_Teams,aux_Score,t+1) + max(0,(lambda1*(2*(aux_Teams-1))/N))*max(0,L_team_rival_team) + lambda1*(aux_sp/N)*max(0,L_team_rival_sp) + lambda2*(aux_sp/N)*max(0,L_team_team_forms);
            end
            
            %%%%%%%%%%F_team
            f_team_own_0 = fun_exp_max(p_sub_f_team_0,q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
            f_team_own_1 = fun_exp_max(p_sub_f_team_1,q_team(aux_Score),L_team(iaux_Teams,min(aux_Score+1,nScore),t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
            
            %computing F^{team,rival team}
            if aux_Teams>1
                aux_F_team_rival_team = @(pl,pf,l) (l/(aux_Teams-1))*(pl*(q_team(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l-1)/(aux_Teams-1))*(pf*(q_team(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
                F_team_rival_team_0 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_0,0);
                F_team_rival_team_1 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_1,1);
            else
                F_team_rival_team_0=0;
                F_team_rival_team_1=0;
            end
            
            %computing F^{team,rival sp}
            if aux_Teams<N/2
                aux_F_team_rival_sp= @(pl,pf,l) ((1-l)/(aux_sp))*(pl*(q_sp(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-1+l)/(aux_sp))*(pf*(q_sp(aux_Score)*F_team(iaux_Teams,min(aux_Score+1,nScore),t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
                F_team_rival_sp_0 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,0);
                F_team_rival_sp_1 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,1);
            else
                F_team_rival_sp_0=0;
                F_team_rival_sp_1=0;
            end
            
            if aux_Teams<N/2
                F_team_team_forms_0 = p_team_f_sp_0*F_team(iaux_Teams+1,aux_Score,t+1,1+0) + (1-p_team_f_sp_0)*F_team(iaux_Teams,aux_Score,t+1,1+0);
                F_team_team_forms_1 = p_team_f_sp_1*F_team(iaux_Teams+1,aux_Score,t+1,1+1) + (1-p_team_f_sp_1)*F_team(iaux_Teams,aux_Score,t+1,1+1);
            else
                F_team_team_forms_0=0;
                F_team_team_forms_1=0;
            end
            
            F_team(iaux_Teams,aux_Score,t,1+0) = ((2*lambda1)/N)*f_team_own_0 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+0) + max(0,(lambda1*(2*(aux_Teams-1))/N))*max(0,F_team_rival_team_0) + max(0,lambda1*((aux_sp)/N))*max(0,F_team_rival_sp_0) + lambda2*((aux_sp)/N)*max(0,F_team_team_forms_0);
            F_team(iaux_Teams,aux_Score,t,1+1) = ((2*lambda1)/N)*f_team_own_1 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+1) + max(0,(lambda1*(2*(aux_Teams-1))/N))*max(0,F_team_rival_team_1) + max(0,lambda1*((aux_sp)/N))*max(0,F_team_rival_sp_1) + lambda2*((aux_sp)/N)*max(0,F_team_team_forms_1);
            
            
            %%%% Storing CCPs
            Psub_f_team(iaux_Teams,aux_Score,t,1)=p_sub_f_team_0;
            Psub_f_team(iaux_Teams,aux_Score,t,2)=p_sub_f_team_1;
            Psub_f_sp(iaux_Teams,aux_Score,t,1)=p_sub_f_sp_0;
            Psub_f_sp(iaux_Teams,aux_Score,t,2)=p_sub_f_sp_1;
            Psub_l_team(iaux_Teams,aux_Score,t)=p_sub_l_team_1;
            Psub_l_sp(iaux_Teams,aux_Score,t)=p_sub_l_sp_0;
            Pteam(iaux_Teams,aux_Score,t,1)=p_team_f_sp_0;
            Pteam(iaux_Teams,aux_Score,t,2)=p_team_f_sp_1;
                   
        end
    end
[k T];    
end

%%
nSim=2500;
Out=zeros(nSim,4);
Teams = zeros(T,nSim);
for k=1:nSim
    state=zeros(T+1,3);
    subs=zeros(T,1);
    mergers=zeros(T,1);
              %[s,n^teams,\ell] 
    state(1,:)=[1,0,0];
    for t=1:T
        rng(t+T*k)
        aux_sub=unifrnd(0,1)<=lambda1;
        aux_Teams = state(t,2);
        aux_sp = N - 2*aux_Teams;
        iaux_Teams=aux_Teams+1;
        aux_Score = state(t,1);
        aux_ell = state(t,3);
        if aux_sub==1
            aux_turn_sp=unifrnd(0,1)<=aux_sp/N;
            if aux_turn_sp==1
                if aux_ell==0
                    aux_turn_leader = unifrnd(0,1)<=1/aux_sp;
                    if aux_turn_leader==1
                        subs(t,1) = unifrnd(0,1)<=Psub_l_sp(iaux_Teams,aux_Score,t);
                    else
                        subs(t,1) = unifrnd(0,1)<=Psub_f_sp(iaux_Teams,aux_Score,t,aux_ell+1);                        
                    end     
                else
                    subs(t,1) = unifrnd(0,1)<=Psub_f_sp(iaux_Teams,aux_Score,t,aux_ell+1);
                end
                if subs(t,1)==1 && unifrnd(0,1)<=q_sp(aux_Score)
                    state(t+1,1) = min(aux_Score+1,nScore);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = 0;
                else
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = state(t,3);                    
                end
            else
                if aux_ell==0
                    subs(t,1) = unifrnd(0,1)<=Psub_f_team(iaux_Teams,aux_Score,t,aux_ell+1);
                else
                    aux_turn_leader = unifrnd(0,1)<=1/aux_Teams;                    
                    if aux_turn_leader==1
                        subs(t,1) = unifrnd(0,1)<=Psub_l_team(iaux_Teams,aux_Score,t);
                    else
                        subs(t,1) = unifrnd(0,1)<=Psub_f_team(iaux_Teams,aux_Score,t,aux_ell+1);                        
                    end                       
                end
                if subs(t,1)==1 && unifrnd(0,1)<=q_team(aux_Score)
                    state(t+1,1) = min(aux_Score+1,nScore);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = 1;
                else
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = state(t,3);                    
                end
            end
        else
            if aux_Teams==N/2
                state(t+1,1) = state(t,1);
                state(t+1,2) = state(t,2);
                state(t+1,3) = state(t,3);
            else
                aux_turn_sp=unifrnd(0,1)<=(aux_sp-(1-aux_ell))/N;
                if aux_turn_sp==1
                    mergers(t,1) = unifrnd(0,1)<=Pteam(iaux_Teams,aux_Score,t);
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2)+mergers(t,1);
                    state(t+1,3) = state(t,3);
                else
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = state(t,3);
                end
            end
            
        end
    end
    Teams(:,k)=mergers;
    %Out(k,:)=[sum(subs),Svalues(state(T,1)),sum(mergers),(sum(mergers)-mergers_data).^2];
end

writematrix(mean(Teams,2),'figure_out_teamstime.csv') 
